Systematically Accelerated Convergence of Path Integrals 
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We present a new analytical method that systematically improves the convergence of path integrals 
of a generic AT- fold discretized theory. Using it we calculate the effective actions S^ p ' for p < 9 
which lead to the same continuum amplitudes as the starting action, but that converge to that 
continuum limit as 1/N P . We checked this derived speedup in convergence by performing Monte 
Carlo simulations on several different models. 
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Since their inception Q path integrals have pre- 
sented an extremely compact and rich formalism for deal- 
ing with quantum theories. They have grown into pow- 
erful tools for dealing with symmetries (including gauge 
symmetry), for deriving non-perturbative results (e.g. 
solitons, instantons, symmetry breaking) and for showing 
connections between different theories or different sectors 
of the same theory (e.g. bosonization, duality) 0, 0- 
They have also consistently allowed us to extend and 
generalize quantization procedures to ever more compli- 
cated systems. Path integrals have brought about a rich 
cross fertilization of ideas between high energy and con- 
densed matter physics by delineating and strengthening 
similarities that exist between statistical and quantum 
systems 0, 0] . Today, path integrals are used both ana- 
lytically and numerically 0, 0, El El El in many other 
areas of physics, chemistry and materials science. They 
are starting to play more prominent roles in several areas 
of mathematics and in modern finance An exten- 

sive review of path integrals and their applications can 
be found in (l3| . 

Unfortunately, we still have very little knowledge of the 
mathematical properties of path integrals. In addition, 
a very small number of path integrals can be solved ex- 
actly. Although the functional formalism has been used 
to derive many general approximation techniques (e.g. 
perturbation, semi-classical expansion, variational meth- 
ods) along with many model-specific approximation tech- 
niques, there still remains a wealth of interesting models 
that can't be analyzed analytically and need to be treated 
numerically. The definition of path integrals as a limit 
of multiple integrals makes their numerical evaluation 
quite natural. The most all-around applicable numerical 
method for such calculations is based on Monte Carlo 
simulations. However, numerical integration of path in- 
tegrals is notoriously demanding of computing time - so 
much so that specific path integral calculations serve as 
benchmarks for new generations of supercomputers. 

Several research groups have focused on improving the 
convergence of path integrals. The best available result 
for a generic theory (valid only for partition functions and 
not for general amplitudes) is the convergence of iV-fold 
discretized expressions as l/N 4 0,EJE3- Related in- 



vestigations have focused either on improvements in short 
time propagation 0, 0, El or have presented model 
specific improvements of the action |2(J, [2l|. Ref. 
gives a useful comparison of several different approaches. 

In order to further significantly speed up numerical 
procedures for calculating path integrals for a generic the- 
ory it is necessary to add new analytical input. In this 
letter we present a systematic investigation of the relation 
between different discretizations of a given theory. A re- 
sult of this investigation is a procedure for constructing a 
series of effective actions having the same continuum 
limit as the starting action S, but which approach that 
limit as 1/N P . We use this procedure to obtain explicit 
expressions for these effective actions for p < 9. 

In the functional formalism the quantum mechanical 
amplitude A(a,b;T) = (b\e~ TH \a) is given in terms of a 
path integral which is simply the N — > oo limit of the 
(N — 1)- fold integral expression 



A N (a,b;T) 



2tt€n 



dqi ■ ■ -dqN-i e 



(1) 



The euclidean time interval [0, T] has been subdivided 
into N equal time steps of length tN — T/N , with qg = a 
and qjy = b. Sn is the naively discretized action of the 
theory. We focus on actions of the form 
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whose naive discretization is simply 
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where S n = q n+ i - q n , and q n = \{q n +i + Qn)- We use 
units in which h and particle mass equal 1. 

The definition of the functional integral makes it nec- 
essary to make the transition from the continuum to 
the discretized theory, a process that is far from unique. 
For theories described by eq. J5J we have the freedom 
to choose any point in [q n ,q n +i] in which to evaluate 
the potential without changing physics - the discretized 
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amplitudes do differ, but they tend to the same contin- 
uum limit. The calculations we present turn out to be 
simplest in the mid-point prescription where V is eval- 
uated at q n . A more important freedom related to our 
choice of discretized action has to do with the possibil- 
ity of introducing additional terms that explicitly vanish 
in the continuum limit. Actions with such additional 
terms will be called effective. For example, the term 
J2 n =o e N $n s(<7n), where g is regular when — > 0, does 
not change the continuum physics since it goes over into 
£ n Jo dt q 2 g(q), i.e. it vanishes as t 2 N . Such terms do 
not change the physics, but they do affect the speed of 
convergence. A systematic study of the relation between 
different discretizations of the same path integral will al- 
low us to explicitly construct a scries of effective actions 
with progressively faster convergence to the continuum. 

We start by studying the relation between 2iV-fold 
and A^-fold discretizations. From cq. Q we see that we 
can write the 2Af-fold amplitude as an A^-fold amplitude 
given in terms of a new action Sn determined by 

— 

e- SN = (^—J 2 J d Xl ---dx N e- S2N , (4) 

where S 2 n is the 2JV-fold discretization of the starting 
action. We have written the 2A r -fold discretized coor- 
dinates Qq, Qi, . . . , Q2N hi terms of q's and x's in the 



following way: Q 2 k = qk and Q 2 k-i = Xk- Note that 
we have go = a , Qn — b, while the N — 1 remaining q's 
play the role of the dynamical coordinates in the A^-fold 
discretized theory. The a;'s are the N remaining interme- 
diate points that we integrate over in eq. |4"|. 

In order to iterate this discretization halving process 
the new action Sn must belong to the same class of ex- 
pressions as the starting actions Sn- It is not difficult to 
show that the naively discretized action does not satisfy 
this requirement, i.e. integration of eq. (jJJ yields new 
types of terms in Sn- In fact, the class of actions closed 
to transformation J3J is of the form 

N-l , S 2 

S N=Y.( 7T- + tNV(q n ) + e N 5 2 l g 1 {q n ) + 

+ e^v $t 92{q n ) + e N $1 g^{q n ) + ■ ■ ■ V (5) 

All the functions appearing above also depend on the 
time step eat. We will not display this dependence ex- 
plicitly in order to have a more compact notation. All 
of the functions are regular in the e^v — > limit making 
these effective actions equivalent to the starting action. 
From eq. Q and JSJ we can derive the following integral 
relation which determines the new action Sn in terms of 
the starting action: 
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The above integral equation can be solved for the sim- 
ple cases of a free particle and a harmonic oscillator, and 
gives the well known results. Note however that for a 
general case the integral in eq. l Efl) i s in a form that is 
ideal for an asymptotic expansion |23l| . The time step cat 
is playing the role of small parameter (in complete paral- 
lel to the role h plays in standard semi-classical, or loop, 
expansion). As is usual, the above asymptotic expansion 
is carried through by first Taylor expanding F(q + y) 
around q and then by doing the remaining Gaussian in- 
tegrals. Assuming that e^r < 1 (i.e. N > T) we have 



All that remains is to calculate the F {2m \q n ys using 
eq. and to expand the potential and all the functions 
gk around the mid-point q n . This is a straight forward 
though tedious calculation. In this letter we illustrate the 
general procedure for calculating S by explicitly giving 
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its expansion to order e^: 
V = V + e N ' 
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In the above relations we expanded V up to e 2 N , gi up 
to ejv, etc. We also disregarded all the higher gk's. The 
reason for this is that the short time propagation of any 
theory satisfies 5„ oc e^v while the term enters the ac- 
tion multiplied by 5^ . In general, if we expand the new 



action S to e^, we need to evaluate only V (up to e^T 1 ) 
and the first p — 1 functions gt (up to e^ -1 k )- The task 
of calculating S to large powers of ej\r is time-consuming 
and is best done with the help of a standard package for 
algebraic calculations such as Mathematica. Using Math- 
cmatica we determined the corresponding expressions for 
p < 9. The above solution for S (and its analogues for 
higher values of p) has the following important property: 
up to O(e^), a coarser TV-fold discretization using S does 
the same job as the 2iV-fold discretization of the starting 
theory i.e. A N (a, b; T) = A 2N {a, b; T) + 0(e%). 

We next iterate the discretization halving process out- 
lined above in order to connect up the 2 s TV-fold and N- 
fold discretizations. Ultimately we will focus on the solu- 
tion when s — > oo, i.e. the one connecting the continuum 
theory with its TV-fold discretization. For the p = 3 case 
that we continue to use as an illustration of the general 
procedure, the above iterative process is governed by: 
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where k = 0, 1, 2, . . . , s — 1. The zeroth iterate corre- 
sponds to the starting action, the last iterate to the ef- 
fective action whose AT-fold discretization is equivalent to 
the 2 s iV-fold discretization of the starting theory. The 
£7v/2 s_fc ~ 1 terms represent the time step of the fc-th it- 
erate in the discretization halving procedure. 

Although the above system of recursive relations is 
non-linear, it is in fact quite easy to solve if we remem- 
ber that the system itself was derived via an expansion 
in ejv ■ Having this in mind we first write all the functions 
as expansions in powers of e^v that are appropriate to the 
level p we are working at. For p = 3, we have 

^^hiB k ' ^ 



V k 

(91) k 

(92) k 



s-fc-1 



F k 



s-fe-l 



El 



(11) 



Putting this into eq. (|10|) we find that Ak+i = Ak- Since 
Aq = V it follows that Ak = V for all k. The remain- 
ing terms obey a linear system of equations that is easily 
solved for given initial conditions. Rather than solving 
the system and then taking the s — + 00 limit of the so- 
lutions, it is easier to set k = s — 1 in the system and 
then take the limit s — ► 00. By doing this we directly 



obtain the continuum limit solution of the discretization 
halving process connecting the continuum theory to its 
TV-fold discretization. The continuum limit solution of 
the p = 3 level system (the solution of the continuum 
limit of the recursive relations given in eq. describ- 
ing discretization halving at the p = 3 level) is: 
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Note that the continuum solution depends only on the 
initial potential V, i.e. it is not sensitive to initial values 
of the gk's as these terms all vanish in the continuum 
limit. In this way we have obtained the effective action 
that gives the best TV-fold discretization of the starting 
theory at the p — 3 level - the one differing from the con- 
tinuum limit by a term of order 1/TV 3 . One can similarly 
obtain continuum limit effective actions for higher 
p. These effective actions lead to TV-fold amplitudes that 
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deviate from the continuum amplitude as 0(e p N ) 

A^(a,b;T)^A( ai b;T) + 0(e p N ). (13) 

It is important to note that one solves for the contin- 
uum limit of the level p system of recursive relations but 
once for all theories, i.e. once the solution is found it 
works for all potentials V. The only requirement for the 
level p solution is that the starting potential is differen- 
tiable 2p — 2 times. Solutions for larger values of p are 
a bit more cumbersome, however, they are just as easy 
to use in simulations. We have found that the growth in 
complexity of the effective actions with increasing p has 
little effect on computation time for p < 4, while simu- 
lations with p = 9 are roughly ten times slower due to 
this. However, this is an extremely small price to pay for 
a gain of nine orders of magnitude in the speed of con- 
vergence. Expressions for effective actions up to p = 9 
can be found on our web site 1241 . 
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FIG. 1: (Color online) Deviations from the continuum limit 
\A^ — A\ as functions of N for p = 1, 2, 4 and 6 for an an- 
harmonic oscillator with quartic coupling A = 10, time of 
propagation T — 1 from a = to b = 1. Nmc was 9.2 ■ 10 9 
for p = 1,2, 9.2 • 10 10 for p = 4, and 3.68 ■ 10 11 for p = 6. 
Dashed lines correspond to appropriate 1/N polynomial fits 
to the data. Solid lines give the leading 1/N behavior. The 
level p curve has a 1/N P leading behavior. 

The analytical derivations presented work equally well 
in both the Euclidean and Minkowski formalism (with 
appropriate ie regularization) , i.e. they are directly ap- 
plicable to quantum systems as well as to statistical ones. 
However, the Monte Carlo simulations used to numeri- 
cally document our analytical results necessarily needed 
to be done in the Euclidean formalism. We analyzed 
in detail several models: the anharmonic oscillator with 
quartic coupling V = | q 2 + |, q 4 and a particle moving 
in a modified Poschl- Teller potential over a wide range 
of parameters. In all cases we found agreement with 
eq. lfT3*]l. Fig. n illustrates this behavior in the case of an 
anharmonic oscillator. We see that the p level data in- 
deed differs from the continuum amplitudes as a polyno- 



mial starting with 1/N P . The deviations from the contin- 
uum limit \A]!p — A\ become exceedingly small for larger 
values of p making it necessary to use ever larger values 
of Nmc so that the MC statistical error does not mask 
these extremely small deviations. For p = 6 we see that 
although we used an extremely large number of MC sam- 
ples (Nmc = 3.68 • 10 11 ) the statistical errors become of 
the same order as the deviations already at N > 8. For 
p = 9 this is the case even for N — 2, i.e. we already get 
the continuum limit within a MC error of around 10~ 8 . 

To conclude, we have presented an algorithm that leads 
to significant speedup of numerical procedures for calcu- 
lating path integrals. The increase in speed results from 
new analytical input that comes from a systematic inves- 
tigation of the relation between discretizations of differ- 
ent coarseness. We have presented an explicit procedure 
for obtaining a set of effective actions that have the 
same continuum limit as the starting action S, but which 
approach that limit ever faster. Amplitudes calculated 
using the iV-point discretized effective action sffi satisfy 
A<fi\a,b;T) = A{a,b;T) + 0(1/N P ). We have obtained 
and analyzed the effective actions for p < 9 and have doc- 
umented the speedup up to 1/N 9 by conducting Monte 
Carlo simulations of several different models. Extension 
to d > 1 is in progress. The derivation of the analogue 
of integral eq. 10 does not seem to present a problem. 
The asymptotic expansion used to solve it is also directly 
generalizable. However, the algebraic recursive relations 
that determine S^ p ' will be more complex and may prac- 
tically limit us to smaller values of p. 
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